/*
 * KineticSolver.h
 *
 *  Created on: 28 Jul 2011
 *      Author: Allan
 */

#ifndef KINETICSOLVER_H_
#define KINETICSOLVER_H_

// C++ includes
#include <vector>
using namespace std;

// GeoReact includes
#include "Multiphase.h"
#include "ReactionSystem.h"
#include "CanonicalReactionSystem.h"
#include "OdeSolver.h"
#include "EquilibriumSolver/TotalEquilibriumSolver.h"

// GeoReact forward declarations
class KineticProblem;

class KineticSolver
{
public:
	KineticSolver();
	
	KineticSolver(const Multiphase& multiphase, const CanonicalReactionSystem& equilibriumReactions, const ReactionSystem& kineticReactions);
	
	~KineticSolver();
	
	void Initialiaze(const Multiphase& multiphase, const CanonicalReactionSystem& equilibriumReactions, const ReactionSystem& kineticReactions);
	
	void SetControlSolution(double epsAbs, double epsRel);
	
	bool Iterate(double& t, double tEnd, double& dt, double T, double P, VectorXd& n);
	
	bool Solve(double& t, double tEnd, double& dt, double T, double P, VectorXd& n);
	
private:
	/// The multiphase object responsible to compute the activities
	Multiphase multiphase;
	
	/// The equilibrium-controlled reactions
	CanonicalReactionSystem eReactions;

	/// The kinetically-controlled reactions
	ReactionSystem kReactions;
	
	vector<Index> jIndexes, eIndexes, kIndexes;
	
	unsigned Nj, Ne, Nk;
	
	VectorXd nj, ne, nk;
	
	VectorXd uj, ujnk;
	
	OdeSolver odeSolver;
	
	TotalEquilibriumSolver eSolver;
	
	/// --- Friend classes --- ///
	friend class KineticProblem;
};

class KineticProblem : public OdeProblem
{
public:
	KineticProblem(const Multiphase& multiphase, const CanonicalReactionSystem& eReactions, const ReactionSystem& kReactions);
	
	void
	SetParameters(double T, double P, const VectorXd& n);
	
	virtual void 
	Function(double t, const VectorXd& ujnk, VectorXd& R);
	
private:
	/// The multiphase object responsible to compute the activities
	Multiphase multiphase;
	
	/// The equilibrium-controlled reactions
	CanonicalReactionSystem eReactions;

	/// The kinetically-controlled reactions
	ReactionSystem kReactions;
	
	/// The temperature (T) and pressure (P)
	double T, P;
	
	/// The composition of the species
	VectorXd n;
	
	/// The number of primary (Nj) and kinetic (Nk) species
	unsigned Nj, Nk;
	
	/// The indexes of the kinetic species (kIndexes)
	vector<Index> kIndexes;
	
	/// The stoichiometric matrices (ukj) and (ukk)
	MatrixXd ukj, ukk;
	
	/// The total equilibrium solver
	TotalEquilibriumSolver eSolver;
};

#endif /* KINETICSOLVER_H_ */
